On  Maximum  Likelihood  Estimators 
of  Shape  and  Scale  Parameters  and  Their 
Application  in  Constructing  Confidence  Contours 


Dl-82-0972 


ON  MAXIMUM  LIKELIHOOD  ESTIMATORS 
OF  SHAPE  AND  SCALE  PARAMETERS  AND  THEIR 
APPLICATION  IN  CONSTRUCTING  CONFIDENCE  CONTOURS 

by 

Sam  C.  Saunders 


Mathematical  and  Information  Sciences  Report  No.  5 
Mathematical  and  Information  Sciences  Laboratory 
BOEING  SCIENTIFIC  RESEARCH  LABORATORIES 


May  1970 


SUMMARY 


If  a  family  of  distributions  is  of  the  form  F[  (t/g)a]  for 
t  >  0,  where  F  is  known  and  a  and  g  are  the  shape  and  scale 
parameters  respectively,  the  maximum  likelihood  estimates  a,  8, 
calculated  from  a  complete  sample  or  an  incomplete  sample  censored 
under  certain  conditions  or  terminated  at  a  given  order  observation, 
have  the  property  that  U  -  a/a  and  V  ■  (g/g)01,  have  distributions 
which  are  parameter-free.  Thus  T  -  U  in  V  has  a  parameter-free 
distribution.  Consequently,  a  tabulation  of  these  distributions 
for  a  specific  sampling  situation  could  be  used  to  perform  the 
usual  statistical  tests  and  to  construct  interval  estimates  for  the 
parameters  in  the  same  way  the  sample  variance,  s,  and  "Student's" 
t-distribution  do  for  the  normal  law. 

Moreover,  it  is  shown  that  this  result  can  be  used  to  obtain 
confidence  contours  along  the  entire  distribution  function  analogous 
to  the  Kolmogorov-Smirnov  bounds  on  the  empiric  cumulative  distribution. 

Taking  the  Weibull  distribution  and  a  sample  of  the  first  3  order 
observations  out  of  5  as  an  example  (in  which  case  the  parameter-free 
properties  of  U  and  V  are  not  new),  computational  procedures  are 
specified  for  determining  the  distribution  of  U  and  V  by  Monte 
Carlo  methods  using  the  latest  random  number  generators.  These 
computing  times  are  given,  as  it  is  for  the  calculations  necessary 
to  compute  a  confidence  contour  along  the  entire  distribution.  Graphs 
of  the  distributions  of  U  and  V  and  the  confidence  contour  for  the 
distribution  are  presented  for  this  case. 


1.  Introduction 


In  a  paper  by  Thoman,  Baine  and  Antle  [12],  It  was  shown  that 
certain  pivotal  functions  of  the  maximum  likelihood  estimators  for 
parameters  in  the  Weibull  family  have  distributions  which  are  parameter- 
free.  (This  fact  was  also  mentioned  in  the  mathematical  appendix  of  the 
proprietary  report  [2].)  'By  using  Monte  Carlo  methods,  this  basic 
result  made  possible  the  production  of  tables  of  the  percentile  points 
of  these  distributions  so  that  confidence  intervals  for  the  parameters 
can  be  determined.  In  turn  this  makes  possible  tests  of  hypotheses 
regarding  the  parameters,  see  [11]. 

It  is  the  purpose  of  this  note  to  generalize  this  result  in 
several  respects.  Firstly,  to  point  out  that  for  any  given  continuous 
distribution  with  support  on  the  positive  line  and  unknown  shape  and 
scale  parameters  the  same  pivotal  functions  of  the  maximum  likelihood 
estimators  are  parameter-free.  Secondly,  that  the  same  conclusion 
holds  for  incomplete  samples  censored  under  certain  conditions  or  termin¬ 
ated  at  a  given  ordered  observation.  Moreover  in  these  cases  it  is 
possible  to  obtain  confidence  contours  along  the  entire  distribution 
function  not  just  of  tall  probabilities  as  given  in  [6].  These  con¬ 
tours  are  not  defective  near  the  tails  as  are  the  Kolmogorov^Smirnov 
bounds  and  by  being  more  specific  are  narrower. 
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2.  The  Model  and  the  Estimates 

Let  R  denote  a  given  survival  distribution  (unity  minus  the 
distribution)  with  support  on  the  positive  real  line.  We  assume  the 
density  f  -  -R'  exists  and  is  differentiable  and  we  denote  the 
hazard  rate  by  q  ■  f/R.  We  now  make  the  assumption 

1°:  The  observable  random  variable  X  has  an  unknown 
survival  distribution  H  within  the  two  parameter 
family  defined  for  given  R  by 

H(x)  -  P[X  >  x)  -  R[(x/B)“]  for  x  >  0, 

where  a  >  0  is  the  shape  parameter  and  B  >  0 

is  the  scale  parameter,  often  called  the  characteristic 

life. 

Alternatively,  we  could  formulate  a  model  with  unknown  scale 
and  location  parameters  and  a  fixed  distribution  by  considering  £nX 
as  an  observable  variate.  Thus  there  is  an  obvious  analogue  for  every¬ 
thing  which  follows  within  that  formulation. 

We  shall  now  introduce  a  sampling  situation  which  often  arises 
in  life  testing,  namely  we  observe  either  the  time  at  which  failure 
takes  place  or  the  time  at  which  the  life  test  is  terminated  with  the 
component  unfailed  and  we  know  which  of  these  events  occurred.  For 
example,  in  fatigue  life  testing  we  are  told  the  total  time  the 
specimen  has  been  fatigued  and  whether  or  not  it  has  broken.  There 


are  similar  situations  in  clinical  trials  in  medicine. 
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Let  X  denote  the  life  length  and  Z  the  random  time  at  which 
the  test  is  terminated  for  any  other  reason  than  failure.  Suppose  we 
observe  the  random  couple  min(X,Z),  £x  <  Z^  where  is  the  indi¬ 

cator  of  the  event  A,  taking  the  value  one  or  zero,  and  min(X,Z)  is 
the  time  the  test  was  operating. 

Let  X^  for  i  »  l,...,n  be  mutually  independent  and  identically 
distributed  life  length  random  variables  with  distribution  F,  density 
F'  and  survival  distribution  F.  Let  Z ^  for  i  »  l,...,n  be  a  set 
of  non-negative  censoring  random  variables  which  may  be  dependent  upon 
the  X's. 

Let  Y±  -  min(Xi,Zi),  I±  -  £xf  <  for  i  -  l,...,n  and  call 

the  random  set  of  Indices  of  failed  items 

A  ■  {i  •  1, . . . ,n:I^  «0}  . 

We  now  make  assumption 

2*:  (Zi>...,Zn)  are  independent,  for  given  A  -  X  of 

{X. :i  {  X}  on  the  event  [x.  >  z  ] . 

i$X  1  1 

This  assumption  means  that  if  a  given  value  of  (x^ . x^)  should 

result  in  a  set  of  (z^ . z^)  then  for  any  1  such  that  >  z^ 

(i.e.,  jeX)  any  larger  observed  value  of  x^  would  have  resulted  in 

the  same  set  of  (z, ,...,z  )  in  distribution.  Thus  for  each  censored 

1  n 

test  the  component  being  even  longer  lived  would  not  alter  the  result 
any. 


•  •  •  •  • 


Lemma:  If  2®  holds  the  density  of  for  i  ■  1 

proportional  to 


n  is 
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(2.1)  TT  F'  (y.)  IT  F(y  ) 
ieX  1  j^X  2 

where  the  constant  of  proportionality  depends  upon  the  failure  observations 
{y^:ieX}  but  not  upon  F. 

Proof:  Consider  the  event 


yl . Yn  " 


V 


-  I 


k 


1, 


k+1 


or  equivalently 

k  n 

(2.2)  [X,  <  Z.][X  -  x,]  [X,  >  Z.  ]  [Z.  -  z.] 

i-1  1  11  1  i-k+1 


I 

n 


-  0] 


Let  g  denote  the  joint  density  of  Z.,...,Z  given  X^,...,X 

n 

assuming  Z. , . . .  ,Z  independent  of  X.  . , . . .  ,X  on  [X  >  Z  ] . 

x  n  K  i  x  ti  ■  •  i  •  x  x 

i-k+1 

The  likelihood  of  the  event  specified  in  (2.2)  is 


Except  for  the  altered  notation  in  (2.1)  this  is  the  result.  || 

Consider  the  special  case  when  the  Z^'s  are  independent  of  the 

X1's.  This  would  include  censoring  at  a  specific  time  determined  by 

either  chance  or  design.  The  case  when  the  Z^'s  are  functionally 

dependent  on  the  X^'s  would  include  the  case  of  censoring  at  a  given 

failure.  Clearly  2°  is  true  in  the  case  of  Z^'s  Independent  of  X^'s. 

th 

To  check  that  2°  is  true  in  the  case  Z^  -  X^  i  ■  l,...,n,  the  k 
ordered  observation,  we  examine  the  joint  density  of  X^j  ,  XJt+1,...,Xn 
on  the  event  <  X±  for  i  -  k+l,...,n.  It  is 
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k(!D  ^‘Sk)’  Jk  f'<v  x(k)<xi 


i  -  k  +  l,...,n 


Note  this  is  the  condition  specified. 

From  (1.1)  we  see  the  density  of  X  is 


(2.3)  f(xa/ga)axa  ^£  a  for  x  >  0  . 


Following  the  notation  introduced  in  the  lemma  we  let  (y^,...,yk) 
denote  the  set  of  observations  of  failed  items,  i.e.,  observations  of 
X  and  (yk+j_»  •  •  •  »yn)  denote  the  set  of  censored  tests,  i.e., 
observations  of  Z. 

From  (2.1),  setting  F' (x)  equal  to  (2.3),  we  may  write  the 
log- likelihood  as 

k  n 

L  “  £  Unf (y“/6a)  +  £n(a/£)  +  (a-l)Jln(y  /£)]  +  £  JlnR(y“/Ba) 

i-1  1  1  i-k+1  1 

Hence 


Z  tn(yi/6)  -  £  (y1/6)V(y1/B)*  k(y“/Ba) 
i-1  i-1 


where  in  general 


<l»l  k(x)  -  q(x)  “  §i  ^  k^q’(x)/q(x)  i,k  -  l,...,n, 


and  similarly 


8L  ka  ,  a  v  /  N a .  .  a/0aN 

3B  -'T  +  3  £  (yi/S)  *i,k<yi/6  >  * 


Thus  the  joint  maximum  likelihood  estimators  a,  3  are  the 
simultaneous  solutions  to  the  equations 


(2.4) 


£  £  (t)  *i,k1(yi/5)“ 


)  •  1 


sWe***5—* 
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(2.5) 


t( j)  ‘»<y1rt)*1,kitf1/S)si  -  Z 


k 

• 

a 


Note  that  if  we  assume  a  Weibull  model,  then  q(x)  ■  1  for 
x  >  0  and  we  have  for  all  i,k  -  l,...,n 

(2.6)  ^  k(x)  -  1  for  all  x  >  0  . 

Naturally  enough  this  results  in  considerable  simplification.  In 

a 

particular  (2.4)  can  be  explicitly  solved  for  6. 

For  the  log-normal  model,  namely 


q(x) 


91 1  (inx) 
x91  (-inx) 


for  x  >  0  , 


where  91  is  the  standard  normal  distribution,  we  find 


(2.7)  i^i>k(x) 


1  +  Enx 
x 

q(x) 


i  S  k 
i  >  k  . 


Substituting  (2.6)  and  (2.7)  into  the  equations  (2.4)  and  (2.5)  we 
find  they  reduce  to  those  given  by  Cohen  in  [4]  for  the  Weibull  and 
log-normal  cases,  respectively. 

x 

Another  case  is  an  extreme  value  distribution  with  q(x)  ■  e  for 
x  >  0.  Here  we  have 


i  >  k 
i  £  k  . 

This  also  results  in  some  simplification  in  (2.4)  and  (2.5). 

Let  W  have  survival  distribution  R  and  for  i  -  l,...,n  set 

(2.8)  wA  »  (y^/6)°  u  -  a/a  v  -  (6/6)°  . 
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Thus  we  see  for  i  ■  l,...,n 

(y1/B)a  -  (wi/v)u,  (y^B)0  -  wi/v 
and  simplifying  (2.4)  and  (2.5)  we  obtain  two  equations  in  (u,v): 

1  n 

k  E  ^^i.k^i5  “ 1 

(2.9) 

k  jE  ^iAn^itk(Ci)  “  k  E^  "  1 

where 

^i  •  (w±/v)u  i  -  l,...,n  . 

It  follows  from  (2.8)  that  W^,...,W^  are  distributed  independently 
of  a  and  8.  We  shall  assume 

3°  (\+i»  •  •  •  »wn)  have  a  ^nown  distribution  independent 

of  o  and  6  where  k  is  the  given  number  of  failed 
(uncensored)  items. 

In  the  case  of  random  sampling  for  a  fixed  time,  it  is  clear 
that  2°  would  not  be  satisfied;  on  the  other  hand,  it  is  clear  that 
sampling  until  a  fixed  number  of  failures  would  satisfy  it.  There  are 
other  situations  which  arise  in  practice  for  which  2°  holds.  For  example, 
in  certain  fatigue  tests  of  structural  details,  failure  takes  place  other 
than  within  the  area  of  primary  concern.  Such  fatigue  may  be  un¬ 
representative  of  service  if  it  is  caused  by  the  abnormalities  of 
local  stress  induced  from  clamping  the  detail  in  the  fatigue  machine. 

The  assumption  is  made  that  the  shape  parameters  of  the  two  distributions 
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of  fatigue  at  the  different  areas  are  the  same  while  the  scale 
parameters  (characteristic  lives)  have  a  known  ratio  determined 
from  the  maximum  deflection  and  the  gross  area  stresses. 

There  follows  immediately  from  equations  (2.9)  the 

Theorem:  If  1°,  2°  and  3®  hold,  then  the  random  variables 
U  *  oi/a  ,  V  ■  (8/8)°  have  a  joint  distribution 
independent  of  a  and  8. 

Also  we  have  the 

Corollary:  If  1°,  2“  and  3°  hold  then  the  random  variable 
T  -  UfcnV  *  a(i,ng  -  Jing)  has  a  distribution 
independent  of  a  and  8. 

The  preceding  theorem  is  a  generalization  of  the  results  given 
in  [12]  for  the  Weibull  distribution.  It  was  not  known  when  the 
survey  article  [7]  was  written  although  it  may  have  been  suspected 
from  sampling  results  for  this  case,  see  the  references  in  [7].  In 
[12]  the  percentiles  of  the  marginal  distributions  of  U  and  V  are 
tabulated  only  for  complete  samples,  for  obvious  reasons.  It  is 
mentioned  there  that  the  parameter- free  property  of  the  pivotal 
functions  holds  when  the  observations  are  censored  at  the  k^  failure. 

Obviously,  if  we  have  a  log-normal  distribution  so  that  (2.7) 
holds  we  obtain  the  theorems  concerning  the  log-normal  model.  These 
can  be  related  to  results  about  the  t-dlstrlbution  through  an  exponential 
transformation.  Thus  the  corollary  above  says  essentially  that  "Student’s" 
result  [10]  about  exact  inference  can  be  extended  to  any  two-dimensional 
family  known  except  for  scale  and  location  parameters  when  using  the 
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maximum  likelihood  estimates. 

The  importance  of  the  preceding  results  is  that  it  establishes 
the  usefulness  of  the  maximum  likelihood  estimators ,  whose  optimum 
properties  are  well  known,  under  a  wide  category  of  censoring  and 
truncation  of  the  data.  Further  it  gives  conditions  under  which  the 
pivotal  functions  of  the  maximum  likelihood  estimators  are  parameter- 
free  and  thus  it  is  possible,  by  varying  R,  to  perform  studies  of 
the  robustness  of  these  estimators. 

If  in  a  given  instance  one  wishes  to  determine  a  joint  con¬ 
fidence  region  in  both  parameters,  or  in  one  parameter  separately  with 
the  other  unknown  (or  the  corresponding  tests)  then  we  need  to  calculate 
the  required  percentile  points  of  the  distribution  which  can  be  done 
using  some  acceptable  Monte  Carlo  procedure.  This  has  already  been 
accomplished  for  certain  cases  for  the  Weibull  distribution  in  [11]. 


3.  Confidence  Contours  on  the  Distribution 


Although  the  need  of  calculating  confidence  contours  along  the 
entire  distribution  function  arises  often  in  reliability  practice 
there  are  but  a  few  methods  of  doing  so. 

Let  H  be  an  estimate  of  the  survival  distribution  Ht3f. 
Previously  in  [9],  we  have  called  ?f  ample  for  H  iff  hTh  *(p) 
for  each  pe(0,l)  has  a  distribution  Independent  of  Hc5f.  (Here 
and  in  what  follows  juxtaposition  of  functions  indicates  composition.) 
For  such  estimators  the  analogue  of  the  Kolmogorov- Smirnov  statistic 

D  -  ^  sup  |H(x)  -  H(x)|  ■  sup  |HH  *(p)  -  p| 
x  p 

and  the  Cramlr-von  Mlses  statistic 

-  1 

-  -n  J  |hT  -  H|2dH  -  n  J  |‘Bh"1(p)  -  p|2dp 

0 

are  distribution-free  with  respect  to  X.  For  the  model  presented 
in  §2  we  obtain 

(3.1)  HH_1(p)  -  R{[VR_1(p)]U)  for  0  <  p  <  1  . 

The  Importance  of  the  probability  integral  transform  to  obtain 
estimates  which  are  parameter-free  is  evident  and  is  not  new;  for 
example,  it  was  studied  in  [5].  Thus  we  see  H(x)  -  R[(x/B)a]  is 
ample  for  any  two  parameter  family  with  R  specified,  since  clearly 
the  distribution  of  the  quantity  HH  *(p)  for  each  pc  (0,1)  does  not 
depend  upon  the  parameters  a, 8. 

It  is  possible  to  obtain  the  distribution  of  D  by  simulation 
and  if  tables  of  the  percentage  points  were  provided  they  could  be 
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used  to  determine  confidence  contours  similar  to  the  well-known 
Kolmogorov-Smirnov  bounds,  e.g.,  see  [3].  However,  we  do  not  favor 
such  bounds  because  of  their  very  bad  behavior  in  the  tails  of  the 
distribution.  This  is  usually  the  region  where  we  are  most  interested 
in  good  behavior.  Instead  we  turn  to  another  method. 

Suppose  that  there  exists  for  each  e<(0,l)  a  continuous 
monotone  Increasing  function,  say  A£,  mapping  (0,1)  onto  (0,1) 
such  that 

(3.2)  P[A~^H  >  H]  -  P[HH_1  >  A£]  -  e 

where  an  inequality  between  functions  indicates  the  inequality  holds 
for  their  functional  values  at  all  points  in  their  common  domain. 

From  (3.2)  an  upper  confidence  contour  of  level  c  for  H  then 
would  be  A~^H.  At  a  later  time  we  mention  the  alterations  necessary 
to  obtain  either  lower  or  two-sided  confidence  contours. 

If  we  substitute  (3.1)  into  (3.2)  and  apply  R-1,  which  is  an 
order-reversing  transformation,  then  take  logarithms  we  obtain 

(3.3)  P[T  +  U  lnR_1(p)  £  inR*1*  (p)  for  all  p*(0,l)]  -  c  . 

If  we  set 

(3.4)  nc(x)  “  *A£R(eX) 

and  make  the  change  of  variable  x  ■  £nR  1  (p>  in  (3.3) 
we  obtain 


c  -  P(T  <  inf(nc(x)  -  xu>]  . 
x 


(3.5) 
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He  now  define  a  functional  on  the  extended  real  line  for  any 
real  valued  function  n  and  any  u  >  0  by 

(3.6)  a)  (u)  ■  inf  [n(x)  -  xu]  . 

n 

— oo<X<°° 

The  question  is,  can  we  find  a  function  n  for  which  u  (U) 

e  ne 

exists  finitely  with  probability  one  and  satisfies  equation  (3.5), 
to  wit, 

(3.7)  e  -  P[T  £  u)  (U)]  . 

'e 

If  so,  it  would  be  theoretically  possible  to  obtain  an  upper 
confidence  contour  along  the  entire  distribution  H  from  knowledge 
of  the  joint  distribution  of  U  and  T.  We  now  seek  conditions  which 
help  with  the  determination  of  ne* 

Let  e  be  the  set  of  Increasing  convex  functions  mapping  the  real 
line  onto  itself  each  of  which  has  a  continuous  derivative  with  the 
positive  real  line  as  its  range.  A  proof  is  easily  given  for  the 

Remark;  If  ntO  then  u  (U)  exists  and  is  a  random  variable 

n 

given  by 

(3.8)  u>n(U)  -  n [ (n * )_1  (U) ]  -  UOi')'1®)  . 

Note  that  if  4>tG  then  mG  where 

(3.9)  n(x)  •  9(x/b)  +  a  for  -  <*>  <  x  <  ® 

for  any  real  a  and  b  >  0  and 

a)  (u)  ■  a  +  u.  (bu)  for  u  >  0. 
n  9 


(3.10) 


13 

It  is  clear  that  by  proper  choice  of  a,b  the  right  hand 
side  of  (3.10)  can  be  determined  for  given  $<<?  so  that  (3.7)  is 
satisfied,  in  theory,  for  each  e  in  the  unit  interval. 

From  the  definition  of  ne  in  terms  of  A£  in  (3.4)  we  find, 
setting  i  *(x)  -  ex 

(3.11)  A  -  RJt-1n  HR”1  . 

£  E 

Thus  the  upper  confidence  contour  of  level  e  is 

(3.12)  A_1H(t)  -  Rr1n“1[oa(t/8)]  for  t  >  0. 

e  e 

In  the  preceding  discussion  an  arbitrary  choice  of  $  was 
made  as  an  illustration  of  the  feasibility  of  the  computation  which 
must  be  done.  Of  course  a  better  choice  of  the  parameterization  could 
be  made  by  considering  the  power  for  the  test  corresponding  to  the 
confidence  contour  for  a  specified  alternative.  It  is  possible  that 
tests  of  an  optimal  nature  could  be  constructed  in  certain  instances 
against  certain  alternatives. 
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4.  Lower  and  Two-Sided  Confidence  Contours 

If  we  are  concerned  with  obtaining  lower  confidence  contours  for 
H  we  would  seek  an  appropriate  monotone  increasing  function,  say  B£, 
for  c  near  unity,  mapping  (0,1)  onto  (0,1)  such  that 

P[H  <  B£H]  -  P [ b”1^  <  H]  -  e. 

Proceeding  as  before  and  making  the  same  change  of  variable  we 
obtain,  in  analogy  with  (3.5),  an  expression  equivalent  with  the 
above,  namely 

(4.1)  P[T  >  sup(v#(x)  -  xU)]  -  E 

E 

X 

where 

(4.2)  v#(x)  -  inR-1B  R(eX)  for  -»  <  x  <  »  . 

E  £ 

Define  for  any  u  >  0  and  any  function  v  for  which  it  exists 

(4.3)  Py(u)  “  sup[v(x)  -  xu], 

x 

If  we  set 

(4.4)  v*(x)  ■  -v(-x)  for  -«  <  x  <  00 

we  see 

(4.5)  ov#(u)  ■  -inflv(x)  -  xu]  ■  -u>v(u)  . 

x 

Thus  if  v«<?  then  v*  is  a  monotone  increasing  concave  function 
mapping  the  real  line  onto  itself  with  a  continuous  derivative  having 
the  positive  real  line  as  its  range.  We  label  this  set  of  functions  C* . 


As  a  consequence  of  (4.5)  we  may  relate  some  of  our  previous  results 
to  this  case.  For  Instance  fix  and  n *G  as  before  In  (3.18)  then 

(4.6)  (n*)  *(y)  -  -n  1(-y)  —  —  bq>  1(a-y) 
and  n  ^*G^. 

If  we  determine  for  a  prescribed  $  a  pair  (a,b)  which  by 
(3.9)  defines  an  upper  confidence  contour  of  level  e,  near  aero, 
as  given  in  (3.7)  in  terms  of  n  then  by  the  use  of  the  trans¬ 
formation  (4.4)  and  the  identity  (4.5)  we  obtain 

(4.7)  P[T  >  a)  (U)  1  -  P[-T  £  o  #(U)J  •  1  -  e. 

-  n  n 

Comparing  (4.1),  we  see,  in  those  situations  where  T  has  a 
distribution  symmetric  about  zero  (e.g.,  when  we  have  complete  samples 
and  the  underlying  survival  distribution  H  possesses  the  requisite 
symmetry)  that  r/  as  given  in  (4.6)  determines  a  lower  confidence 
contour  of  level  1  -  e  which  is  near  unity. 

In  the  more  usual  case  where  we  do  not  have  -T  with  the  same 
distribution  as  T  a  determination  must  be  made  of  the  appropriate 
(a,b)  in  (4.6)  in  order  to  satisfy  (4.1).  A  discussion  of  the 
feasibility  of  this  we  defer  until  later. 

To  obtain  both  upper  and  lower  confidence  contours  for  H 
simultaneously  we  seek  A  and  B  both  monotone  increasing  functions 
mapping  the  unit  interval  onto  itself  such  that 

e  -  P[BH  i  H  2  AH]  -  P[A-1H  >  H  i  B_1H]  . 

From  the  first  equality  above,  we  obtain  as  before 
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e  -  P[iR-1B(p)  <  T  +  UtR-1Ap,  for  all  p«(0,l)] 

(4.8)  -  P[pv#(v)  £  T  *  a)n(v)] 

where  la  defined  In  terms  of  B  by  (4.2),  n tC  Is  defined 

in  terms  of  A  by  (3.11)  and  p.w  were  defined  respectively  by 
(4.3)  and  (3.6). 

From  (4.4)  we  seek  both  v,n *6  for  which 

(4.9)  c  -  P[-o>v(U)  <  T  <  con (U) ]  . 

The  determination  of  an  appropriate  v  and  n  in  (4.9),  or  the 
corresponding  n  in  (3.7)  or  in  (4.1)  can  be  done  rather 

straightforwardly  if  we  fix  and  choose  the  appropriate  values 

of  the  parameters  (a,b)  across  a  two  dimensional  subspace  of  6  as 
defined  in  (3.9). 

If  ve  do  so  in  the  two  sided  case  we  have  from  (3.10) 

e  -  P[-w.(b*U)  -  a#  £  T  i  a  +  uA(bU)] 

with  the  obvious  use  of  notation.  If  we  exercise  two  degrees  of  freedom 
and  set  b*  -  b,  a*  ■  a  we  obtain 

(4.10)  e  ■  P ( | T |  <  a  +  w.  (bU)] 

from  which  the  proper  choice  of  the  parameters  may  be  more  easily  made. 

A  method  of  choosing  the  appropriate  parameters  a,b  we  take  up  next. 


5.  Monte  Carlo  Techniques 


This  discussion  Is  to  make  the  point  that  with  the  latest 
random  number  generators  available  and  the  computing  capability 
currently  extant,  the  distribution  of  the  relevant  statistics  from 
U,V  can  in  practice  be  determined  by  simulation  with  as  much  precision 
and  speed  as  could  be  done  in  many  instances  using  numerical  calculation 
were  the  distributions  available  in  closed  form.  Moreover,  the  wide 
variety  of  sampling  situations  which  can  arise  in  practice  makes  it 
virtually  impossible  to  provide  tables  of  more  than  limited  usefulness 
for  even  one  specific  choice  of  R. 

He  now  outline  a  method  for  the  determination  of  either  a  con¬ 
fidence  interval  on  a  with  6  unknown  (or  B  with  a  unknown  or 
both)  or  confidence  contours  along  the  entire  distribution  (or  on  one 
side  only)  for  the  Welbull  failure  model  when  a  given  number  of  ordered 
observations  have  been  obtained.  The  general  procedure  may  be  inferred 
from  this  particular  one. 

Let  y^,...,y^  be  the  first  ordered  observations  from  n  i  k 
independent  machine  geuerated  exponential  variates  with  unit  mean. 
Recall,  for  the  Welbull  model,  R(x)  »  e  x.  Solve  for  u  in  the 
equation  x(u)  "0  where 


and 


by  using  Newton-Raphson  iteration  procedure. 
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Now  compute  for  determined  as  a  solution  of  X(u)  -  0  the 

value 


The  value  is  one  observation  of  (U,V).  This  process 

can  be  repeated  m  times  from  ndependent  samples  of  y's  to  obtain 
an  independent  sample  of  size  m  from  (U,V).  The  value  of  m  can 
be  made  sufficiently  large  to  determine  the  joint  distribution  of  U,V 
(or  the  marginals)  to  a  degree  of  accuracy  limited  by  the  machine  pro¬ 
cedure  that  generates  random  numbers.  The  desired  percentiles  of  this 
empiric  distribution  are  then  tabulated,  which  can  then  be  used  to 
obtain  confidence  intervals  or  regions  in  the  obvious  way.  We  do  not 
pursue  this  matter  farther  because  of  the  similarity  with  work  done 
previously  in  [6].  But  as  an  example  we  present  the  empirical  marginal 
distributions  of  U  and  V  in  Figures  1  and  2,  respectively,  for  k  -  3, 
n  ■  5  with  m  -  4500.  Time  for  both  computations  was  5.67  minutes  on  the 
IBM  360. 

The  generated  random  numbers  recommended  here,  and  used  in  [2], 

are  of  the  type  called  composite  congruential  generators.  These  second 

generation  methods  appear  to  be  better,  that  is  they  satisfy  more 

stringent  statistical  tests  of  randomness,  than  those  of  the  simple 

congruential  generators  used  previously.  In  the  particular  method 

utilized,  three  generators  are  mixed  for  the  IBM  360  each  of  which  will 

32 

produce  a  full  period  of  residues  relatively  prime  to  the  modulus  2  . 

30 

Consequently,  these  mixed  generators  will  produce  2  distinct  random 
numbers  before  repeating.  This  method  is  presented  in  detail  and  its 
practical  advantage  discussed  in  [8].  To  obtain  our  exponential 
observations  with  unit  mean  we  merely  take  the  negative  of  the  natural 
logarithm  of  the  uniformly  distributed  observations  generated  by  the 


Figure  1.  Distribution  of  U  ■  a/a,  where  a  is  the  shape  parameter 

of  the  Weibull  distribution  and  a  is  the  maximum  likelihood 
estimate  based  on  the  first  three  observations  failing  out  of 
five. 
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mixed  congruentlal  method. 

We  new  outline  the  construction  of  confidence  contours  for  the 
Weibull  distribution  by  supposing  <J>  was  fixed  and  finding  the 
appropriate  a,b  from  the  two-dimensional  subset  of  C  so  that 
either  (3.7),  (4.7)  or  (4.9)  is  satisfied.  Recall  that  T  ■  OinV 
and  let 


r(a,b)  ■  P[T  i  a  +  ai  (bU)  ] 

or  we  might  replace  T  by  |t|  as  in  (4.10). 

A  feasible  procedure  then  is  to  take  a  mesh  of  values  in  the 

half  plane  a,b  >  0  and  on  the  basis  of  our  sample  (U, ,V. ,V  ) 

m  m 

calculated  from  the  appropriate  random  number  generators,  compute  the 
relative  frequencies  of  toe  occurrence  of  the  appropriate  events.  We 
tabulate 

r(vV  *  -  &  & s  *i +  Wk'3 

where  as  before  Is  the  indicator  of  the  relation  it  being  one  if 

true  and  zero  otherwise.  From  the  evidence  of  these  results  we  may 
wish  to  interpolate  to  find  more  appropriate  values  and/or,  as  we 
proceed,  to  increase  the  sample  size  m  and/or  refine  the  mesh. 

As  a  test  case  we  chose 


4>(x) 


from  which 


-in(l-x) 


for  x  >  0 
for  x  £  0 


Vu> 


u  -  1  -  u£nu 


if  u  >  1 
if  u  <  1 


1  -  u  +  £nu 
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A  mesh,  presented  in  Table  1,  for  the  range  of  values 

a  -  -.6(.l)-.3  and  .4(.1).8 
b  -  . 25 ( . 25) 2 . 75 


with  k  ■  3,  n  ■  5  and  m  ■  4,500  was  computed  on  the  IBM  360  in  5.09 
minutes.  (It  was  repeated  to  verify  that  all  entries  were  accurate  to 
two  significant  ‘igures.)  As  an  example  the  entry  a  ■  .8,  b  ■  .25 


has  r(a,b)  -  .9484  thus 


n  1(y) 


bln(l  +  y  -  a) 

b(l  -  ea-y> 


a  s  y  <  » 

<  y  s  a 


can  be  used  in  conjunction  with  (3.12)  to  obtain  an  upper  confidence 
contour  of  level  .95.  A  graph  is  presented  in  Figure  3. 
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15 


Figure  3.  Upper  confidence  contour  of  level  .95  for  a  Weibull  distribution 
with  a  ■  2,  3  ■  20  using  a  ■  1.776,  6  ■  21.36  calculated 

from  the  first  three  ordered  observations  out  of  five,  with 
contour  function  $  defined  by  4>(x)  •  ex  -  1  for  x  >_  0, 

<t> (x)  ■  -£n(l-x)  for  x  <  0. 
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